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Abstract 

This paper shows how to find lower bounds on, and sometimes solve globally, a 
large class of nonlinear optimal control problems with impulsive controls using semi- 
definite programming (SDP). This is done by relaxing an optimal control problem 
into a measure differential problem. The manipulation of the measures by their 
moments reduces the problem to a convergent series of standard linear matrix in- 
equality (LMI) relaxations. After providing numerous academic examples, we apply 
the method to the impulsive rendezvous of two orbiting spacecrafts. As the method 
provides lower bounds on the global infimum, global optimality of the solutions can 
be guaranteed numerically by a posteriori simulations, and we can recover simulta- 
neously the optimal impulse time and amplitudes by simple linear algebra. 



1 Introduction 

Optimal control problems are still an active area of current research despite the availability 
of powerful theoretical tools such as Pontryagin's maximum principle or the Hamilton- 
Jacobi-Bellman approach, that both provide conditions for optimality. However, numeri- 
cal methods based on such optimality conditions rely on a certain number of assumptions 
that are often not met in practice. In addition, state constraints are particularly hard to 
handle for most of the methods. 

On the other side, many numerical methods have been developed that deliver locally 
optimal solutions satisfying sufficient optimality conditions. However, the users of these 
methods are often left to wonder if a better solution exists. For example, in the particular 

^CNRS; LAAS; 7 avenue du colonel Roche, F-31077 Toulouse; Prance. 

^Universite de Toulouse; UPS, INSA, INP, ISAE; UTl, UTM, LAAS; F-31077 Toulouse; France 
^Faculty of Electrical Engineering, Czech Technical University in Prague, Technicka 2, CZ-16626 
Prague, Czech Republic 

^Institut de Mathematiques de Toulouse, Universite de Toulouse; UPS; F-31062 Toulouse, Prance. 



case of impulsive controls, it is often not known if more regular solutions could provide 
a better cost. For a recent survey on impulsive control see e.g. [ill ^^"^ references 
therein. See also [7j for a recent application and for more reference^ For historical works 
see e.g. [13 [ISl 122] and also 

This paper presents a method based on [121 121] but covering a larger class of solutions, 
including impulsive controls. This algorithm provides a sequence of non-decreasing lower 
bounds on the global minimizers of affine-in-the-control polynomial optimal control prob- 
lems. In particular, it may assert the global optimality of local solutions found by other 
methods. As importantly, the algorithm is also able to provide numerical certificates of 
infeasibility or unboundedness for ill-posed problems. Finally, in some cases, it is also 
possible to generate the globally optimal control law. 

At the end of the paper, this method is successfully applied to the problem of coplanar 
space fuel-optimal linearised rendezvous. We show with two different examples from the 
literature that the proposed algorithm is able to retrieve the impulsive optimal solution 
conjectured by running a direct approach based on the solution of a Linear Programming 
(LP) problem. Without assuming the nature of the propulsion (continuous or impulsive), 
the obtained impulsive solution is certified to be a global fuel-optimal solution. 

1.1 Contributions 

The paper improves the method presented in [T21 [21] in the following ways: 

• Impulsive control can now be taken into account. 

• Because controls are represented by measures and not by variables, the size of SDP 
blocks is significantly reduced. This allows to handle larger problems in terms of 
number of state variables as well as to reach higher LMI relaxations. 

• Total variation constraints can be handled very easily. 

These three improvements make it altogether possible to tackle problems such as con- 
sumption minimization for space rendezvous, the other significant contribution of this 
paper. 

1.2 Notations 

Integration of a function / : — )■ M with respect to a measure /z on a set X C is 
written /^/(x) d/i(x). The Lebesgue or uniform measure on X is denoted by A whereas 
the Dirac measure concentrated at point x is denoted hj 6x- A measure /x is a probability 
measure whenever J = 1. The support of measure is denoted by supp /i. The 
indicator function of set X (equal to one in X and zero outside) is denoted by Ix- 

^We are grateful to Terence Bayen for pointing out this reference to us. 
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F{X) is the space of Borel measurable functions on X, whereas BV{X) is the space of 
functions of bounded variation on X. M.[z] is the ring of polynomials in the variable z. 
B{X) denotes the Borel a-algebra associated with X. 

If k eN"' denotes a vector of indices then x'^ with a; G M" is the multi-index notation for 
Uxi'. The de gree of the index k is deg k = ^ki. Finally, is the set of all indices for 
which deg k < e N*^. 



2 The optimal control problem 

This paper deals with the following nonlinear optimal control problem 

Vixo) = inf I(xo,u)= [ h(t,x(t))dt+ [ H(t)u(t)dt+ hTixiT)) 
«WeF([o,T])- Jo Jo 

s.t. x{t) = f{t,x{t)) + G{t)u{t), te[0,T] ^'^^ 

x{0) =Xoe Xo, x{T) e Xt, x{t) G X C M". 

where the dot denotes differentiation w.r.t. time and the prime denotes transposition. 
Criterion I{xo, u) is affine in the control u, and V is called the value function. It is assumed 
that all problem data are polynomials, meaning that all functions are in ]R[t, x], and that 
all sets are compact basic semialgebraic. Recall that such sets are those which may be 
written as {x : aj(x) > 0, ■? = !,..., m} for some family {ai}^^^ C ]R[x]. A mild technical 
condition (implying compactness of X) must be satisfied [I3l Assumption 2.1], but it 
is often met in practice (for instance, an additional standard ball constraint ^ < 
enforces the condition). The reason for making these assumptions will be apparent in the 
later sections. 

Without additional assumptions and constraints, the infimum in problem ([T]) is generally 
not attained in the space of measurable functions [26j. For this reason, in this paper we 
consider problems for which controls are allowed to be generalized functions, i.e. measures, 
thereby extending the original formulation as follows: 

Vr(xo) = inf I(xo,w)= [ h(t,x(t))dt+ [ H(t)dw(t) + hTixiT)) 

«,WeBy([o,T])- Jo Jo 

s.t. dxit) = fit, x{t)) dt + G{t)dw{t), t G [0, T] 

x{0) = Xo G Xo, x(T) G Xt, x{t) G X C M" 

(2) 

where Vr stands for the relaxed value function. In particular, in problem ^ controls 
may be impulsive: the (vector) control can been seen as a (vector) distribution of the first 
order and it is therefore the distributional derivative dw{t) of some (vector) function of 
bounded variation w{t) G BV{[0,T])"', see e.g. [22] and ^ §4] or also [51, Prop. 8.3]. 
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3 The measure problem 



In this section, we formulate problem ([2]) into a measure differential problem, a necessary 
step towards obtaining a tractable SDP problem. Optimal control problems involving 
measures have been introduced to accept solutions that are ruled out or ill-defined in 
classical optimal control, see e.g. [26]. Multiple solutions, impulsive or chattering controls 
can be handled naturally by the associated measure problem. This section, rather than 
providing rigorous proofs, outlines the main ideas behind this transformation. 

A few remarks are worth pointing out. First of all, it is crucial that G{t) be a matrix of 
smooth functions, an hypothesis automatically fulfilled by polynomials. As a matter of 
fact, multiplying distributions with such functions is a well-defined operation (unlike e.g. 
the product of two distributions). Therefore, except for some very particular cases [IS], 
G cannot be a function of states Xj that could potentially present jump discontinuities. 
To simplify notations, we have simply assumed that G depends on t onljj^ Secondly, in 
the absence of impulses, the distributional differential is the traditional differential, and 
the dynamics are classical differential equations with controls dw{t) = u{t)dt which are 
absolutely continuous with respect to the Lebesgue measure. Finally, it must be noted 
that state trajectories x{t) are themselves functions of bounded variations, being the sum 
of two such functions, and that this is their broadest class in the sense that there is no 
more general distribution solutions for the states [22j. 

Because distributional derivatives of functions of bounded variation on compact supports 
can be identified with measures [T8| §50], the dynamics in problem ([2]) may be inter- 
preted as a measure differential equation. As X C M*^ is assumed to be compact, by 
one of the Riesz representation theorems [TUl §36.6], these measures can be put in dual- 
ity correspondence with all continuous functions v{t,x{t)) supported on [0,T] x X. We 
will use these test functions to define linear relations between the measures. Note that 
because continuous functions on compact sets can be uniformly approached by polynomi- 
als by virtue of the Stone- Weierstrass theorem, it is enough to consider polynomial test 
functions v{t,x{t)) G M[[0,T] x X]. 

By Lebesgue's decomposition theorem [TUl §33.3], we can split the control measures w{dt) 
into two parts: their absolutely continuous parts with density u : [0,T] — )• (with 
respect to the Lebesgue measure) and their purely singular parts with jump amplitude 
vectors Ut^ G supported at impulsive jump instants tj, j G J, with J a subset of 
Lebesgue measure zero of [0,T], not necessarily countabl^ We write 

w{dt) = u(t)dt + ^G{tj)ut^6t^{dt) 

to model jumps in state-space 

x+{tj) = x-{tj) + G{tj)ut^, Vj G J. 

■^In all rigour, it could be possible to include state jumps in G, but this requires a careful definition of 
what is meant by integration, as done e.g. for studying stochastic differential equations. This goes well 
beyond the scope of this paper. 

^We suspect however that for the control problems studied in this paper, subset J can be assumed 
countable without loss of generality. 
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Now, given an initial state xq G Xq and given a control w{t) G i?\^([0, T])™, denote by 
x{t) G BV{[0,T])"' the corresponding feasible trajectory. Then for smooth test functions 
V : [0,T] X M" ^ M, it holds 

T 

dv{t,x{t)) = v{T,x{T)) -v{0,x{0)) 



\dt \dx J 



Gudt 



dx 

+ X^^(^i>a;+(tj)) -v{tj,x~{tj)). 

We are going to express the above temporal integration ([s]) along the trajectory in terms 
of spatial integration with respect to appropriate and so-called occupation measures. For 
this purpose, define: 

• The time-state occupation measure 

fx[xo,w{t)]{Ax B) = [ lB{xit))dt, \/AeB{[0,T]), ^B e B{X) 

J A 

which measures the occupation of A x i? by the pair {t, x{t)) all along the trajectory. 
Note that we write fi[xo, w(t)] to emphasize the dependence of fi on initial condition 
Xq and control w{t). However, for notational simplicity, we may use the notation 
fi. By a standard result on Borel measures on a cartesian product, the occupation 
measure fi can be disintegrated into 



fi{AxB)= [ 

J A 



dt. 



where ^{dx\t) is the distribution of x G M", conditional on t G [0,T]. It is a 
stochastic kernel, i.e., 

— for every t G [0,T], 1 1) is a probability distribution on X, and 

— for every B G B{X), C,{B | ■) is a Borel measurable function on [0,T]. 

In our case, since the initial state Xq and the control w{t) are given, the stochastic 
kernel ^{dx \ t) is well defined along continuous arcs of the trajectory as 

e(i?|t) = /B(x(t)) = 4w(5), Wte[0,T]\J, WBeB{X). (4) 
On the other hand, at every jump instant tj G J, we let 

" \(lx-(t,),x+(tj)]) ' ' ' 

This means that the state is uniformly distributed along the segment linking the 
state before and after the jump, the above denominator ensuring that |t) has 
unit mass for all t. 
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• The control-state occupation measure 

p[xo,w{t)]{Ax B) = f i{B\t)dw{t), VA e i3([0,T]), \/BeB{X). 

J A 

• The final state occupation measure 

liT[xo,nj{t)]{B) = Ib{x{T)), \JB G B{Xt). 

With these definitions, Eq. ([3]) may be written in terms of measures as: 

v{T,x) duT^x) — v{^,Xq) = (5) 



1 + (I)'/) ir^oMt,^)- 

[0,T]xX V^^ \OX J J J[0,T]xX \OX J 

Lt + ^")) ^^^^ I ^) dt + J2v{t„x^{t,))-v{t„x-it,)). 

Similarly, the criterion in ^ to evaluate the trajectory and the control reads: 
I{n,iy,HT)= / hdfx+ H du + hrdnx- 

J[0,T]xX J[0,T]xX Jxt 

In view of the above formulation with occupation measures, one may now define a relaxed 
version (or weaA; formulation) of the initial (measure) control problem ([2]). First note that 

Vr{xo) = inf I{fx[xo,w{t)],u[xo,w{t)],n[xo,w{t)]) 

w(t) 

where the infimum is taken over all the occupation measures defined above, corresponding 
to a given initial condition xq and control w{t). Second, instead of searching for a control 
w{t), we search for a triplet of measures that solves the infinite dimensional problem: 

VMixo) = inf I{ij[xo],u[xo\,ht[xo\) 

under the trajectory constraints ^ for all f G ]R[t, x] and the support constraints supp fx = 
supp z/ = [0, T] X X, supp /iT = Xt- The measures now depend only on initial condition 
xq, since they just have to satisfy linear constraints This motivates the notation n[xo\, 
h'lxo], i^tIxo] in the above problem. This problem is an obvious relaxation of problem ^ 
which is itself a relaxation of ([T]), hence 

Vm{xo) < Vnixo) < V{xo). 

In the remainder of the paper, we will deal with this relaxed version of the occupation 
measures problem. However, for a well-defined control problem (|2| one expects that in fact 
Vm{xo) = Vr(xo) and that an optimal solution of the relaxed problem will be the triplet 
of occupation measures corresponding to an optimal trajectory of problem ^ with given 
initial state Xq and control w{t). Note that for the standard polynomial optimal control 
problem ([T|, without impulsive controls, and under additional convexity assumptions, it 
has been proved in [12j that indeed Vm{xo) = Vr{xo) = V{xq). 
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3.1 Initial state with a given distribution 

Recall that the occupation measures defined in the previous section all depend on xq. 
Observe that if /io is a given probability measure on Xq C M" and if one now defines: 

fi{A X B) = J^^ n[xo]{A X B) rf/io(xo), 
u{A X B) = J^^ iy[xo]{A x B) rf/io(xo), 
Ht{B) = /^^ /iT [2:0] (5) d/io (xo) 

for all A G B{[0,T]) and B G B{X), then 



/(/i[;Uo],z/[/io],/iT[H) = / I{l^[xo],iy[xo],HT[xo])diJ^o{xo) 

JXo 

becomes the expected average cost associated with the trajectories and with respect to 
the probability measure /iq on Xq. 

Therefore, the relaxed problem with measures now reads as follows: 
VA/(/io) = inf /(/i[/io], z^[/io]>/^T[H) = / hdj2+ / Hdu + / hrd/j^T 

M,!^,/iT J J J 

s.t. [ vdjjiT— I V dj^Q = [ ( TT + 1 o~ ) f] djj, + [ (7:— ) Gdv 



dt \dx J J J \dx J 

supp yU = supp V = [0, T] X X, supp fiT = Xt- 

Note that in this case, the stochastic kernel ^{dx\t) along continuous arcs of the trajectory 
is generally not a Dirac measure as in Q, unless yUo is a Dirac measure supported at xq 
and the optimal control w is unique. 

By solving this relaxed problem we expect that its optimal value satisfies 



Vm(/^o) = / VM{xo)dij,o{xo), 

JXo 

i.e. that Vm(/^o) is the expected average cost associated with optimal trajectories and 
with respect to the probability measure on Xq. 



3.2 Free initial state 

In this case, in addition to the control we also have the freedom of choosing the best 
possible initial state. For this purpose introduce an unknown probability measure /zq on 
Xq. Then the relaxed problem with measures, analogue of ([6]), reads almost the same 
except that: 

• we now optimize over /x, z/, /io, /^t; 

• in the support constraints we introduce the additional constraint supp ^0 = Xq. 

By solving this relaxed problem we now expect that its optimal value denoted Vm{Xq) 
satisfies: 

VAf(Xo) = inf VAf(/io) = inf Vm{xq). 
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3.3 Decomposition of control measures 



All measures in ^ are positive measures, except for the signed measures u which deserve 
special treatment for our purposes. Using the Jordan decomposition theorem ^U, §34], 
these measures may be split into a positive part and negative part z/~, that is z/ = 
— z/~, both being positive measures. 

This decomposition has the added benefit of providing an easy expression for the Li 
norm of the control, which is sometimes to be constrained or optimized in some problems. 
Indeed, define the total variation control measure by 

lul = u'^ + . 

The total variation norm of the measure u is just the mass of i.e.. 



WWtv = j d\v\ 



3.4 Handling discrete control sets 

It is often desirable to restrict the set of admissible controls to be a subset of M. Here we 
will limit ourselves to the very important case of handling discrete control sets. Let us 
assume that controls u are only allowed to take their values in f/ = {ui, Um}- Define Vi 
as the probability measures of choosing controls Ui. Then clearly, the total probability of 
choosing one of the controls in U must be 1 at each time along the trajectory. Then the 
control measures v are simply the linear combination of the probability measures weighted 
by their respective control values: 



i 

Using the same method as in ^ we have Wv{t) G ]R[t]: 

Note that with this substitution, all measures involved in the measure problem are now 



positive; there is no need to implement the trick of §3.3[ Using these extra constraints, it 
is now possible to solve bang-bang control problems. 



3.5 Summary 

To summarize, the advantages for introducing the relaxed control problem (|6| with mea- 
sures are the following: 

• controls are allowed to be measures with absolutely continuous components and 
singular components including impulses; 
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state constraints are easily handled via support constraints; 

the initial state has a fixed given distribution on some pre-specified domain; 

a free initial state in some pre-specified domain is also allowed. 



4 The moment problem 



So far, the hypothesis of polynomial data has not been used, but its crucial importance 
will appear in this section, where measures will be manipulated through their moments. 
This will lead to a semi-definite programming (SDP) problem featuring countably many 
equations. 

Define the moments of measure /z as 



y^,= / z'dii{z). (7) 
Jx 

Then, with a sequence y = (yk), k G N", let Ly : M.[z] — )■ M be the linear functional 



k 



Define the moment matrix of order (i G N associated with y as the real symmetric matrix 
Md (y) whose (z, j)th entry reads 

MMi,j] = Ly = Vz, J G Nl 

Similarly, define the localizing matrix of order d associated with y and h G M.[z] as the 
real symmetric matrix Md{hy) whose {i,j)th entry reads 

Md{hy)[i,j] = Ly {h{z) z^^^) = J2 hky^+J+k, Vz, J G N^. 

k 

As a last definition, a sequence y^ = (y^) is said to have a representing measure if there 
exists a finite Borel measure /i on X, such that relation ([T]) holds for every k G N". 

Now comes the crucial result of the section: a sequence of moments y'^ has a representing 
measure defined on a semi-algebraic set = {x : af(a;) > 0, i = 1, 2, . . .} if and 
only if Mdiyf") h 0, Vrf G N and Md{a'^ y^') ^0, Vrf G N and Vaf defining set Xf" [131 
Theorem 3.8]. This has the very practical implication that the measure problem defined 
in ^ has an equivalent formulation in terms of moments. Indeed, because all problem 
data were assumed to be polynomial, the criterion in ^ can be transformed into a linear 
combination of moments to be minimized: 

= inf m'y^ + {bn'y" + {h^n'y^^ = b'y (8) 

y 
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where the infimum is now over the aggregated sequence y of moments of all the measures. 
Because the test functions were also restricted to be polynomials, the constraints in ^ 
can be turned into countably many linear constraints on the moments: 

A^y^ + A'y'' + A^"'y^"' + A^'^y''^ =Ay = Q. (9) 

The only non-linear part are the SDP constraints for measure representativeness, to be 
satisfied Vrf G N: 

M,(y'^°)^0, M,(af'y^°)^0, (10) 
M,(y'^-) ^0, M,(ary^^) ^0. 



5 LMI relaxations 

The final step to reach a tractable problem is relatively obvious: we simply truncate 
the problem to its first few moments. Let di G N be the smallest integer such that all 
criterion monomials belong to ^2^"^. This is the degree of the so called first relaxation. 
For each relaxation, we reach a standard LMI problem that can be solved numerically by 
off-the-shelf software by simply truncating Eq. (|8]), ^ and (10) to involve only moments 
in N2j'^, with d > di the relaxation order. 

Observe that dj > di ^ V^.^ > V^j . Therefore, by solving the truncated problem for 
ever greater relaxation orders, we will obtain a monotonically non-decreasing sequence of 
lower bounds to the true cost. In the examples below, we will see that in practice, the 
optimal cost is usually reached after a few relaxations. 



6 Academic examples 

In this section, many examples are presented to showcase the different features of the 
method. Ex. [T] to |5] are variations of the same basic problem to give a thorough tour 
of the method's capabilities. Ex. [6] is taken from the literature and shows how the 
method compares to, or rather nicely complements, existing optimal control algorithms. 
All examples use GloptiPoly [9] for building the truncated LMI moment problems and 
SeDuMi 123] for their numerical solution. 



Before proceeding to the examples, define the marginal Md {y, z) of a moment matrix with 
respect to variable z as the moment matrix of the subsequence of moments concerning 
polynomials of z only. 

Example 1 (Basic impulsive problem). 

V = m{[ x\t)dt 
Ht)Jo 
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Figure 1: Trajectory for Ex. [T] 



such that 



x{0) = 1, x{2) 
x^{t) < 1. 



1 
2 



In this introductory example, it is straightforward to notice that the optimal solution 
consists of reaching the turnpike x{t) = by an impulse at initial time t = 0, and 
likewise, departing from it by an impulse at final time t = T = 2, see Fig. [1} 

The associated measure problem reads: 



Vm = inf / X dfi 

1^'" J[0,T]xX 



such that 



dv 
di 



d/i + 



Xt JXo J[0,T]xX 

Ho = 5o Xo = {1} fiT = 5i Xt 
X = {x eR: 1- x^ >0} . 



[0,T]xX 



dv 
dx 



dv ^veR[t,x] 



Using the procedure outlined above, one obtains a series of truncated moment problems 
that can be solved by semi-definite programming. Letting y^- = j fx^ d^, the first LMI 
relaxation is 

Vm = inf Z/02 



11 



subject to the hnear constraints associated to the dynamics: 







= 


y^^ 


-yt^ 


= 2/00 


y'ol 


- y'^l 


= 2/oo ~ Voo 


ll'r 

yw 




= 2<o 


i/ll 


-y'll 


= 2/01 + 2/10 -2/10 


yo2 


- y',i 


= 2^01 — "^yoi , 



to the SDP representativeness constraints for r = {/i, z/^, }: 



2/oo 2/io 2/oi 
2/io 2/20 2/11 
2/01 2/11 2/02 



^ 0, ?/oo - 2/02 > 0, 



and to the boundary conditions: 

2/fo" 2/o^r 2/2^0° 2/rr 2/0^2"] = [1 01001], 



2/ro" 2/o7 2/2^0^ 2/ff 2/0^21 = [1 2 ^ 4 1 |] . 

It turns out that the optimal value Vm = is estimated correctly (to numerical toler- 
ance) from the first relaxation on and that the optimal trajectory x{t) = can easily 
be recovered. Indeed, the marginal Md{y^-,x) is the length of the time interval multi- 
plying a truncated moment matrix of a Dirac measure concentrated at x = 0, while its 
marginal with respect to t equals a truncated Lebesgue moment matrix on the [0, 2] in- 
terval. More importantly, one can recover the optimal controls as the marginal M^iy", t) 
is the weighted sum of Dirac measures located at the impulse times, the weights being 
the impulse amplitudes. In summary, we can recover numerically the optimal measures 

fi{dt,dx) = I[Q^2]{dt)So{dx), h'{dt,dx) = —6o{dt)I[Q^{dx) + 62{dt)I^Qi-^{dx). 

Example 2 (Total variation constraints). We take back Ex. [Tjwith an additional con- 
straint on the total variation of the control: 



/ \u{t)\ dt<l 
Jo 



whose measure equivalent reads: 

II'^IItv < 1- 

Clearly, the solution of Ex. [l| with a total variation of |, violates this extra constraint, 
so the algorithm should converge to another solution. Again, from the first relaxation 
on, the cost of the associated truncated moment problem is |. It is also plain to see 
that Md (|/^,x) is the truncated moment matrix of a Dirac located at x = |, hinting a 
trajectory x(t) = |. On the control side, starting from the second relaxation, it also 
becomes evident that {y" it) is the truncated moment matrix of the signed measure 
~f ^0 + I ^2, revealing the times and amphtudes of impulses compatible with admissible 
controls. This leads to the trajectory of Fig. |2| which therefore is an optimal solution of 
the problem. 
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Figure 2: Trajectory for Ex. [2] 




Figure 3: Trajectory for Ex. [3] 



Example 3 (Discrete control set with chattering). We take again Ex. [T] with the ad- 
ditional constraint that the control takes its value in the set U = {±1}, using the 
method explained in Section |3.4[ The solution to this problem is easy to infer: reach 
the turnpike x(t) = as quickly as possible by applying the negative control until t = 1, 
then chatter with equal probability to remain on the turnpike until ^ = |, after which the 
positive control must be applied until t = 2 (see Fig. |3|. This solution has an optimal 
cost of I ~ 0.375. Compare this value with those of Table [l| which presents the evolution 
of the criterion with respect to the relaxation order of the truncated problem. After the 
fourth relaxation, the marginal w.r.t. x of the control measure corresponding to the con- 
trol u(t) = +1 closely approaches the positive measure |/[]^ |](dx) + /[| 2]('ia;) while the 
marginal w.r.t. x of the control measure of u(t) = —1 converges to /[o,i]((ix) + ^ I^^ a^^dx), 
as expected. 

Example 4 (Infeasible problem). If the problem is infeasible, it may be detected by the 
infeasibility of one of the LMI relaxations. Take Ex. [T] with the additional total variation 
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Table 1: Criterion as a func 


ion of LMI relaxation 


d 


1 


2 


3 


4 


yd 


0.000 


0.288 


0.368 


0.372 



order for Ex. [3] 



constraint 1^(^)1 dt < \ that puts the end point out of reach from the starting point. 
Indeed, at the first relaxation, the LMI problem is flagged as infeasible with a Farkas dual 
vector, providing a certificate of infeasibility of the original problem. 

Example 5 (Unbounded problem). If the problem is unbounded, it will be detected at 
the first LMI relaxation. Consider the problem of maximizing the total variation of a 
linear control problem: 

sup / \u{t)\dt (11) 

u{t) Jo 

such that 

x{t) = u{t) 

a;(0) = 0, x(l) = 

x'^{t) < 1. 

As expected, the LMI problem from the first relaxation on is flagged as unbounded because 
its dual is flagged as infeasible. 

Example 6 (Bang-bang control of the Vanderpol equation). Consider the following time- 
optimal problem of the Vanderpol equation: 

inf T 

u(t)GU 

such that 

Xi{t) = X2(t) 

Mt) = -Mt) - i^lit) - 1) Mt) + «W 

x{0)=[-OA -0.6]', x(T) = [0.6 0.4]' 
U = {±1}. 

In [23], this problem is solved by applying a gradient-based optimization technique on 
a parametrization of admissible trajectories, with a minimum time of 2.14. However, 
this method can only prove the local optimality of solutions. Applying our method, we 
obtain a cost of 2.15 at the fifth relaxation, providing a (numerical) certificate of global 
optimality for that local solution. 



7 The fuel-optimal linear impulsive guidance rendezvous 
problem 

In this section, the proposed approach is applied to the far-range rendezvous in a linearised 
gravitational field. This problem is defined as a fixed-time minimum-fuel impulsive orbital 
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transfer between two known circular orbits. Under Keplerian assumptions and for a 
circular rendezvous, the complete rendezvous problem may be decoupled between the 
out-of-plane rendezvous problem for which an analytical solution may be found [S] and 
the coplanar problem. Therefore, only coplanar circular rendezvous problems based on 
the Hill-Clohessy- Wiltshire equations and associated transition matrix pj are considered 
for numerical illustration of the proposed results. The general framework of the minimum- 
fuel fixed-time coplanar rendezvous problem in a linear setting is recalled in [5] and [1] 
where an indirect method based on primer vector theory is proposed. Considering the 
necessity of easy-to-implement numerical solution for on-board guidance algorithms, direct 
methods based on linear programming (LP) problem may be used as in [16]. For an a 
priori fixed number of impulsive manoeuvres and using a classical transcription method [3] 
[H] , the genuine infinite-dimensional problem may be converted into a finite-dimensional 
approximation given by the following LP problem: 



V, 



LP 



S.t. 



N 

min V'llueJIi 

u ' 

i=l 



N 



x{ei) = xo, x{9f) = Xf 



:i2) 



i=l 



where $ is the Hill- Clohessy- Wiltshire transition matrix, i? = [ 02x2 I2 ]' and ue^ is the 
vector of velocity increments at 6i in the local vertical local horizontal (LVLH) frame [Ij. 
Time has been changed to the true anomaly 6 for the independent variable as is usual in 
the literature [S], and it ranges in the interval [^1, 6f]. Note that this formulation implies 
that only the impulsive solution of the general linear rendezvous problem may be obtained 
for a fixed number of velocity increments. 



To be consistent with our previous notations we let t 
impulsive optimal control problem ^ writes 



6, 6(] 



and e 



f 



T. Our 



V, 



M 



S.t. 



inf 

w{t) 



dx 



x(0) 



\dwi\{t) + \dw2\{t) 



'0 





1 


0" 




'0 


0" 











1 


x{t)dt + 

















2 


1 





_0 


3 


-2 










1 






x{T) 




Xf 







dw{t) 



where state components model positions (X, Z) = {xi,X2) in the orbital plane, and their 
respective velocities {X,Z) = {x3,X4). It clearly encompasses formulation (12) since 



it allows to consider continuous or impulsive thrusters as well. In both cases, the fuel 
consumption is measured by the one-norm of vector function [6*1, 6*/] — )• [20] whereas 

the two- norm of this vector is used in general in the literature, see [3], [T] and references 
therein. 

For the sake of comparison between these two approaches, two academic examples taken 
from [S] are presented. 
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Example 7 (In-plane rendezvous 1). Consider the first case presented in [5]. It consists of 
a coplanar circle-to-circle rendezvous with zero eccentricity. The rendezvous manoeuvre 
must be completed in one orbital period with boundary conditions Xq = [ 1 ]' 
and X/ = [ ]'. This type of rendezvous is usually difficult to handle by numer- 
ical methods because of its singularity due to the high number of symmetries involved. 

With a grid of = 50 points, the LP algorithm gives a two-impulse solution at the initial 
and final times of the rendezvous without interior impulse nor initial or final coasting pe- 
riod. The optimal impulses are both horizontal and opposite Mo = —U2tt = [ 0.05305 ] . 
The fuel cost is given by Vlp = 0.1061. The LMI method has no difficulty to recover 
the optimal solution given by the LP algorithm. A cost of 0.1061 is obtained for each 
relaxation. It is then easy to extract from the matrices that the optimal solution for the 
first control consists of two symmetric impulses of magnitude 0.0531 at the initial and 
final times, while the second control is identically 0. The optimal trajectory in the orbital 
plane is depicted in Figure |4] where -|- indicates the 50 points of discretization. 




-0 25' 1 1 ' 1 1 1 1 

-0.2 0.2 0.4 0.6 0.6 1 1.2 

Figure 4: Trajectory in the orbital plane (X, Z) in LVLH: Case 1 of [3] 
Figure [5] shows position, velocity and impulses history versus true anomaly. 
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anomaly [rad] 




Z 3 4 

anomaly [rad] 



-0.05 




Z 3 4 

anomaly [rad] 



Figure 5: Positions (X solid, Z dashed), velocities {X solid, Z dashed) and impulses (on 
X axis): Case 1 of 



Example 8 (In-plane rendez-vous 2). As a second example, the third case of |5j| is revis- 
ited. The rendezvous is nearly identical to the previous one except for the final condition 
that imposes to reach the target with relative velocity of 0.427 in the Z direction, namely 
xo = [ 1 ]' and a;/ = [ 0.427 ]'. 

Again, a grid of = 50 points is used when running the LP algorithm. It converges to 
a four-impulse trajectory depicted in Fig. |6| The numerical results are summarized in 
Table 1 
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0.3 r 



O.Z - 



0.1 - 



- 



-0.1 - 



-0.2 - 



-0.3 - 




-0 4 I 1 1 1 1 1 1 1 1 

-0.2 0.2 0.4 0.6 0.6 1 1.2 1.4 

x[m] 

Figure 6: Trajectory in the orbital plane (X, Z) in LVLH: Case 3 of [5] 



Table 2: Impulse times and amplitudes for Ex. [8] 



LMI method 


LP method 






M2 











-0.0386 








-0.0392 





1.791 


+0.109 





1.795 


+0.109 





4.495 


-0.109 





4.488 


-0.109 





6.283 


+0.0389 





6.283 


+0.0392 






Using our algorithm, we reached the same criterion (within numerical tolerance) after the 
fourth relaxation (see Tab. [3]). As usual, the controls can be inferred from the moment 
matrix of the z/ measures. Indeed, ui converges to the measure J2i'^di)i^9i with impulse 
amplitudes (weji and anomaly 9i taken from Table [2| while z/2 converges to an all zero 
measure. Not only does this result prove the global optimality of the conjectured solution 
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within the class of all impulsive solutions no matter the number of impulses, but it also 
shows that it is optimal over all measure thrust solutions. 



Table 3: Criterion as a function of LMI relaxation order for Ex. [8] 



d 


1 


2 


3 


4 


Yd 


0.0463 


0.0680 


0.2188 


0.2972 



Finally, position, velocity and impulses history are illustrated in Figure [7} Note the 
symmetry of the optimal four-impulse solution. 




Z 3 4 

anomaly [rad] 
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Figure 7: Positions (X solid, Z dashed), velocities {X solid, Z dashed) and impulses (on 
X axis): Case 3 of |5] 
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8 Conclusion 



The focus of this work is on actual computation of optimal impulsive controls for systems 
described by ordinary differential equations with polynomial dynamics and polynomial 
(semialgebraic) constraints on the state. State trajectory and controls are measures which 
are linearly constrained, resulting in an infinite-dimensional linear programming (LP) 
problem consistent with the formalism of our GloptiPoly software |9j. This LP problem 
on measures can then be solved numerically via a hierarchy of linear matrix inequality 
(LMI) relaxations, for which off-the-shelf semi-definite programming (SDP) solvers can be 
used. The optimal impulse sequence can then be retrieved by simple linear algebra, and 
global optimality can be verified by a posteriori simulation or comparison with suboptimal 
control sequences computed by alternative techniques. 

For space rendezvous, our technique can be readily adapted to cope with state (e.g. 
obstacle avoidance) constraints, as soon as they are basic semialgebraic. Other criteria 
than the total variation can also be handled. Smoother solutions can be expected, maybe 
consisting of a mix of absolutely continuous and singular controls, including impulsive 
controls. 
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